Surperhard monoclinic BC6N allotropes: First-principles investigations
Qu Nian-Rui1, Wang Hong-Chao1, Li Qing1, Li Yi-Ding1, Li Zhi-Ping1, †, Gou Hui-Yang2, Gao Fa-Ming1, ‡
Key Laboratory of Applied Chemistry, Yanshan University, Qinhuangdao 066004, China
Center for High Pressure Science and Technology Advanced Research, Beijing 100094, China

 

† Corresponding author. E-mail: zpli@ysu.edu.cn fmgao@ysu.edu.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 21671168 and 21875205), the Hebei Natural Science Foundation, China (Grant No. B2015203096), and the Qinhuangdao Science and Technology Support Program, China (Grant No. 201703A014).

Abstract

Via structural searching methodology and first-principles calculations, we predicted two new BC6N allotropes, a C-centered monoclinic BC6N (Cm-BC6N) and a primitive-centered monoclinic BC6N (Pm-BC6N). The lattice vibrations, elastic properties, ideal strength, theoretical hardness, and electronic structure of the predicted BC6N were investigated systematically. Our results reveal that Cm-BC6N is more favorable energetically than graphite-like g-BC6N above 20.6 GPa, which is lower than the transition pressures of r-BC6N, t-BC6N, and Pm-BC6N. Both Cm-BC6N and Pm-BC6N are indirect semiconductors with band gaps of 2.66 eV and 0.36 eV, respectively. Cm-BC6N exhibits the excellent ideal shear strength of 53.9 GPa in (011)[01 ], much greater than that of Pm-BC6N (25.0 GPa in (010)[101] shear direction), and Cm-BC6N shows a much lower anisotropy in shear strength than Pm-BC6N. The Vickers hardness of Cm-BC6N is estimated to be above 80 GPa, which is more outstanding than those of t-BC6N and r-BC6N.

PACS: 62.20.-x
1. Introduction

Due to its broad industrial applications, superhard material is always of great scientific interest.[1,2] As we all know, the diamond and cubic boron nitride (c-BN) are well-studied superhard materials with excellent physical properties, such as wide band gaps, high hardness and thermal conductivity, large shear and bulk modules. However, the adhibition of diamond has been conditioned by its own imperfections, including brittleness, reacting with iron easily, and generating oxidizing reaction in the air at a high temperature. Cubic boron nitride shows better performance than diamond above, but its hardness and mechanical properties are obviously weaker than those of diamond.[3] B–C–N ternary system is considered to be harder and stronger than c-BN, chemically and thermally more stable than diamond.[4] Therefore, searching superhard materials from the B–C–N ternary system has been a focus of scientific research. Zhao et al.[5] successfully synthesized stoichiometric c-BC2N and c-BC4N crystals with diamond-like structures by using a diamond-anvil cell. Then Fan et al.[615] did systematic investigations of these structures based on the first-principles calculations, such as crystal structure, band gap, optical spectrum, ideal strength and hardness, in order to establish facts and reach new conclusions.

As we know that B–C–N compounds can consist of the typical sp3 covalent bonds containing B–C, B–N, C–N, and C–C bonds, the mechanical properties thus would be better with the increase of carbon content.[1618] Graphite-like BC6N (g-BC6N)[19] has been found as an n-type semiconductor, which was synthesized by the chemical vapor deposition. The diamond-like t-BC6N and r-BC6N were proposed by Luo et al.,[20] which were obtained by replacing two carbon atoms in the eight-atom cubic diamond unit cell with a nitrogen atom and a boron atom. According to their results, the Vickers hardness of the two BC6N phases is about 80 GPa, higher than that of c-BN (66 GPa) and c-BC2N (76 GPa). In fact, the first-principles and their evolution methods are always used in predicting possible crystal structures or interpreting the mechanical properties of superhard materials.[21,22]

In the present work, a C-centered monoclinic BC6N (Cm-BC6N) and a primitive-centered monoclinic BC6N (Pm-BC6N) are proposed by structure searching method based on first-principles calculations. The structural, electronic, and mechanical properties of the proposed BC6N allotropes are investigated systematically. The Cm-BC6N and Pm-BC6N have dynamical and mechanical stability in ambient circumstance. We not only find their semiconductor nature from the band structures, but also can draw a conclusion that Cm-BC6N is more favorable than the previous predicted t-BC6N and r-BC6N in thermodynamics. By computing the ideal strength and theoretical hardness, we find out the superhard nature structures.

2. Calculation method

Crystal structures of Cm-BC6N and Pm-BC6N were explored by using the evolutionary methodology through particle swarm optimization as implemented in CALYPSO package,[23,24] which has successfully predicted various structures.[2528] Subsequent structural relaxation and property calculations were performed based on the density functional theory[29,30] with the Vienna ab initio simulation software package (VASP) code, with a projector augmented wave (PAW) scheme[31] to show the electrons–ion interaction. The generalized gradient approximation (GGA) with the Perdew–Burke–Ernzerhof (PBE)[32] functional was employed to describe the exchange and correlation potentials. From our text data, a plane-wave kinetic energy cutoff of 500 eV and quantitative value of k-points sampling of 0.25 Å−1 in Brillouin zone were used to give good convergence of the total energies. Optimization of the internal atomic positions and the lattice parameters was achieved by the minimization of forces and stress tensors. After calculating the elastic constants, the bulk modulus and shear modulus were figured out following the Voigt–Reuss–Hill (VRH) approximations.[33] From the local density approximation in plane wave basis, the phonon dispersions were calculated by the density functional perturbation theory, which has been implemented the Quantum ESPRESSO code[34] with the Troullier–Martins (TM) pseudopotentials.[35] Static enthalpy differences calculated relative to the graphite-like g-BC6N structure for the various structures are plotted in Fig. 2 as a function of pressure. The calculations of stress–strain relations were performed following the method of Refs. [3638], which have calculated the strength of several strong solids successfully.[39,40] By the micro hardness models for covalent crystal, we acquired the hardness.[4143]

3. Results and discussion

The optimized structures of Cm-BC6N and Pm-BC6N are shown in Figs. 1(a) and 1(b). Both of Cm-BC6N and Pm-BC6N have irregular six membered rings along the [100] and [010] directions. There are B–N bonds in Cm-BC6N structure, which is the key distinction to Pm-BC6N. All atoms are four coordinated with sp3 hybridization, resulting in six membered rings like chair/boat. More detailed structure information of Cm-BC6N and Pm-BC6N is listed in Table 1. For comparison, the structure parameters of t-BC6N and r-BC6N are also provided. The unit cell of Cm-BC6N contains sixteen atoms. Two B, two N, and four C atoms locate at four inequivalent crystallographic 2a sites, respectively. The other eight C atoms averagely locate at two inequivalent crystallographic 4b sites. The unit cell of Pm-BC6N contains eight atoms. One B and three C (from C1 to C3) atoms locate at four inequivalent crystallographic 1a sites, respectively. One N and the other three C (from C4 to C6) atoms locate at four inequivalent crystallographic 1b sites, respectively. As shown in Table 1, the mass density of Cm-BC6N (3.446 g/cm3) is the largest among the BC6N allotropes, suggesting that the spatial structure of Cm-BC6N is the most compact one among them.

Fig. 1. Structures of (a) Cm-BC6N and (b) Pm-BC6N; the yellow, orange, and green spheres represent B, C, and N atoms, respectively.
Table 1.

The equilibrium lattice parameters of a, b, c, density ρ, volume, and atomic positions.

.

To assess the thermodynamic stability of the BC6N (Cm-BC6N, Pm-BC6N, t-BC6N, and r-BC6N) phases, the relative enthalpies (without temperature effect) are calculated at the pressure range of 0–60 GPa. The pressure–enthalpy relationship calculated for all BC6N phases is shown in Fig. 2. At zero pressure, g-BC6N is the most stable phase, and Cm-BC6N, Pm-BC6N, t-BC6N, and r-BC6N are higher than g-BC4N by 0.42 eV, 0.75 eV, 0.60 eV, and 0.44 eV per atom, respectively. The enthalpy of Cm-BC6N is the lowest among the other phases within the entire pressure range, suggesting that Cm-BC6N should be accessed more easily than Pm-BC6N, t-BC6N, and r-BC6N. Under hydrostatic compression, Cm-BC6N becomes more favorable energetically than g-BC6N above 20.6 GPa; and it has a lower transition pressure than 22.2 GPa for r-BC6N, 35.8 GPa for t-BC6N, and 55.5 GPa for Pm-BC6N. Moreover we calculate the phonon dispersion to prove the dynamical stability of Cm-BC6N and Pm-BC6N. As shown in Fig. 3, there are no imaginary frequencies observed in the whole Brillouin zone, which indicates that the two phases keep dynamically stable at zero pressure.

Fig. 2. Enthalpy as a function of pressure for BC6N allotropes relative to g-BC6N.
Fig. 3. The phonon-dispersion curves of (a) Cm-BC6N and (b)Pm-BC6N at 0 GPa.

By the Hookeʼs law, we calculate the single-crystal zero-pressure elastic constants Cij to comprehend Cm-BC6N and Pm-BC6N in elastic stability and mechanical properties.[44] For monoclinic symmetry, there are thirteen independent single-crystal elastic constants. As shown in Table 2, the calculated elastic stiffness constants Cij of Cm-BC6N and Pm-BC6N satisfy the mechanical stability criteria for a monoclinic structure,[45] indicating the elastically stability at ambient condition.

Table 2.

Calculated zero-pressure elastic constants Cij (GPa), bulk modulus B (GPa), shear modulus G (GPa), Youngʼs modulus E (GPa), and Poissonʼs ratio υ of Cm-BC6N and Pm-BC6N.

.

Superhard materials not only manifest high bulk and shear moduli but also exhibit low Poissonʼs ratio, because of hardness strongly depending on the plastic deformation.[41] Additionally, the Voigt–Reuss–Hill approximation is used to reckon the polycrystalline elastic moduli as the Cij have been calculated, and the results are listed in Table 2. The calculated bulk modulus of Cm-BC6N (395.3 GPa) within GGA in the present work is the greatest one, followed by r-BC6N (393.3 GPa), t-BC6N (383.7 GPa), and Pm-BC6N (373.7 GPa). Note that the bulk moduli of t-BC6N and r-BC6N are 410.4 GPa and 399.9 GPa in Ref. [20]. The order of shear modulus is the same as that of bulk modulus: Cm-BC6N (462.3 GPa) r-BC6N (458.3 GPa) t-BC6N (452.2 GPa) Pm-BC6N (423.0 GPa). When describing the directional degree of covalent bonds about a material, we always choose another momentous parameter called Poissonʼs ratio υ. All of the proposed phases have low Poissonʼs ratio (0.08–0.09), which indicates the possibilities of directional covalent bonding. Furthermore, all the BC6N phases exhibit the brittle behavior due to the low ratio of B/G (0.85–0.88) according to Pughʼs rule.[46]

The Youngʼs module (E) is one of the most important parameters which measures the stiffness of a solid material. We show the three-dimensional (3D) representation and the corresponding two-dimensional (2D) projection of the Youngʼs modulus for Cm-BC6N and Pm-BC6N and in Fig. 4. The shape of the 3D curved surface of Cm-BC6N is close to a sphere and the 2D projections in the , and xy planes are close to a circle, which show that Cm-BC6N has slight anisotropy. The same thing happens to Pm-BC6N. As can be seen in Fig. 4, the , and xy plane areas of Cm-BC6N are larger than those of Pm-BC6N, indicating that the Youngʼs modulus property of Cm-BC6N is superior to that of Pm-BC6N.[47]

Fig. 4. The 3D representations (left panel) and 2D projections (right panel) of Youngʼs moduli for (a) Cm-BC6N and (b) Pm-BC6N. Note that the negative sign only denotes the direction corresponding to the positive one.

The calculated electronic band structure and density of states (DOS) of Cm-BC6N and Pm-BC6N at 0 GPa are presented in Fig. 5. It can be seen from the left panel of Fig. 5 that Cm-BC6N and Pm-BC6N are semiconductor with an indirect band gap of 2.66 eV and 0.36 eV, respectively. The band gap of Cm-BC6N is larger than that of the previous proposed r-BC6N (2.09 eV) and t-BC6N (2.06 eV). The band gap of Pm-BC6N (0.36 eV) is the smallest among them. From the data, Cm-BC6N has the top of the valence band at point and the bottom of the conduct band at Z point. Pm-BC6N has the top of the valence band at point and the bottom of the conduct band at E point.

Fig. 5. Electronic band structures (left panel) and density of states (right panel) of (a) Cm-BC6N and (b) Pm-BC6N. The energy is measured from the valence top.

To further understand the chemical bonding and electronic properties, we calculate the partial density of states (PDOS) of Cm-BC6N and demonstrate the primary valence-electron states of B, C, and N atoms. From the right panel of Fig. 5, it can be find out that there are harmonic 2p–2s overlaps between B(C, N)-2s and B(C, N)-2p states over the whole valence bands in both Cm-BC6N and Pm-BC6N structures, which ascertains the hybridization of the s and p orbitals and the strong interaction between B, C, and N atoms. It is in accordance with the bonding behaviors in the previous reports of B–C–N system. In Cm-BC6N, the B(C, N)-2s mainly contributes to the valence bands below about 10.0 eV, and the B(C, N)-2p mainly contributes to the valence bands between about 10.0 eV and 0 eV. We note the nonmetallic nature of Cm-BC6N by the large energy band gap which is clear and near the Fermi level.

We also calculate the ideal strength to connect the atomic level strength under plastic deformation and permanent large strain.[48,49] The calculated tensile stress versus tensile strain of Cm-BC6N, Pm-BC6N, t-BC6N, and r-BC6N are shown in Fig. 6. As shown in Fig. 6(a), the stress responses of Cm-BC6N in the seven principal directions with the peak tensile stresses are within the range of 55.9–118.6 GPa. The anisotropy ratio of the ideal tensile strengths for Cm-BC6N is . Nevertheless, the weakest tensile strength of Cm-BC6N is 55.9 GPa in the [011] direction at a strain of 0.08. As shown in Fig. 66(b), the strong stress responses in the seven principal directions with the peak tensile stresses are between 32.9 GPa and 91.8 GPa for Pm-BC6N. The anisotropy ratio of the ideal tensile strengths for Pm-BC6N is . Nevertheless, the weakest tensile strength of Pm-BC6N is 32.9 GPa in the [001] direction at a strain of 0.05. As shown in Fig. 6(c), the strong stress responses in the five principal directions with the peak tensile stresses are between 52.3 GPa and 170.6 GPa for t-BC6N. The anisotropy ratio of the ideal tensile strengths for t-BC6N is . Nevertheless, the weakest tensile strength of t-BC6N is 52.3 GPa in the [111] direction at a strain of 0.08. From Fig. 6(d), we find that with the peak tensile stresses, the strong stress responses in the three prime directions are within 82.8–50.9 GPa for r-BC6N. The anisotropy ratio of the ideal tensile strengths for r-BC6N is . Nevertheless, the weakest tensile strength of r-BC6N is 50.9 GPa in the [210] direction at a strain of 0.08. From the data above, the weakest tensile strength of Cm-BC6N is the largest among the BC6N allotropes, whereas the weakest tensile strength of Pm-BC6N is the smallest. Both of the ideal tensile strengths of Cm-BC6N and Pm-BC6N are significantly larger than that of β-BC2N.[50] 7

Fig. 6. Calculated tensile stress versus tensile strain for (a) Cm-BC6N, (b) Pm-BC6N, (c) t-BC6N, and (d) r-BC6N in principal symmetry directions.
Fig. 7. Calculated shear stress versus shear strain for (a) Cm-BC6N and (b) Pm-BC6N in principal symmetry directions.

We also investigate the stress–strain relations of Cm-BC6N and Pm-BC6N in various principal symmetry directions under shear deformation. The anisotropy ratio of the ideal shear strengths for Cm-BC6N is . The weakest shear strength of 53.9 GPa of Cm-BC6N is found along the (011)[01 shear deformation path. The anisotropy ratio of the ideal shear strengths for Pm-BC6N is . The weakest shear strength of 25.0 GPa of Pm-BC6N is found along the (010)[101] shear deformation path. The magnitude of the ideal shear strength of Cm-BC6N approaches to 58.3 GPa found for the (111)[11 slip system of c-BN.[51] Besides, it is found that the lowest peak shear stresses of Cm-BC6N and Pm-BC6N are lower than their lowest peak tensile stress, respectively, indicating that the mechanical failure modes of these two BC6N allotropes are dominated by the shear type.

Due to the connection of the hardness and the bulk modulus or shear modulus, we can acquire the hardness. The hardness calculations are performed based on the microscopic models of hardness constructed by Gao et al.[4143] and Chen et al.[52,53] respectively, which have been widely applied currently. In Gaoʼs microscopic hardness models, the hardness of a covalent material is determined by two factors: the properties of bonding at equilibrium and the number of bonds per unit area. The later includes the bond length, charge density, iconicity, etc. The Vickers hardness of the μ-type bond can be calculated as follows: where is the number of valence electrons of the μ-type bonds per cubic angstroms, and it is determined by the number of effective valence electrons per μ bond and the bond volume , which can be expressed as . is Phillipsʼs homopolar band gap of the μ-type bonds,[54] which characterizes the strength of the covalent bond.[41] is ionicity of chemical bond in crystal scaled by Phillips,[55] and , where and .

In practical applications, we define the geometric average hardness of all kinds of chemical bonds in crystals as follows: where is the number of the μ-type bonds composing the actual complex crystal.

We calculate the atomic bond hardness of BC6N and show the result in Table 3, and the calculation results of c-BC4N, β-BC2N, c-BN, and diamond are also given for comparison. The theoretical Vickers hardness values of r-BC6N (80.2 GPa) and t-BC6N (80.0 GPa) are according with the calculation results of Ref. [20]. The theoretical Vickers hardness values of c-BC4N(80.3 GPa) and β-BC2N(75.5 GPa) agree well with experimental results of Refs. [15] and [6]. This proves the reliability. The hardness values of Cm-BC6N and Pm-BC6N are predicted to be 80.6 GPa and 77.4 GPa by Gaoʼs model and to be 84.0 GPa and 76.5 GPa by Chenʼs model, respectively. Coincidently, the hardness by Chenʼs model is compared favorably with for the two novel BC6N allotropes. Both methods indicate that Cm-BC6N is a superhard material with the high Vickers hardness above 80 GPa, which is evidently larger than that of c-BN (66.4 GPa), β-BC2N (75.5 GPa), and c-BC4N (80.3 GPa); furthermore, Cm-BC6N is harder than r-BC6N (80.2 GPa) and t-BC6N (80.0 GPa); nevertheless, it is slightly inferior to diamond (93.8 GPa). The results show the superior mechanical properties of Cm-BC6N and Pm-BC6N in potential technological applications.

Table 3.

Hardness calculations of Cm-BC6N, Pm-BC6N, t-BC6N, r-BC6N, c-BC4N, c-BN, and diamond, where is the bond length, is the electron density expressed in number of valence electrons per cubic angstroms, is the homopolar gap, is the valence-conduction band gap, is the ionicity of bond, is geometric average of all bonds vickers hardness, compared to the hardness of (Chen) and the experimental hardness , respectively.

.
4. Conclusion

In summary, two surperhard BC6N with monoclinic structures were proposed. The structural, electronic, and mechanical properties were investigated by first-principles calculations. Cm-BC6N and Pm-BC6N are both mechanically and dynamically stable, and Cm-BC6N is suggested to be more easily synthesized than Pm-BC6N, t-BC6N, and r-BC6N. We found that the lowest stress peak of Cm-BC6N is 53.9 GPa, which is higher than that of other general hard materials. Comparing with other BCN allotropes which have been proposed and Pm-BC6N, Cm-BC6N has a higher hardness from Gaoʼs model (80.6 GPa). In Cm-BC6N, the B and N atoms connect each other by the prominent and strong C–C bond which results in its superior mechanical properties. Both of Cm-BC6N and Pm-BC6N are significantly harder than β-BC2N and c-BN, suggesting the intrinsic superhard materials.

Reference
[1] Haines J Leger J Bocquillon G 2001 Annu. Rev. Mater. Res. 31 1
[2] Liu A Y Cohen M L 1989 Science 245 841
[3] Wentorf R H Jr. 1957 J. Chem. Phys. 26 956
[4] Horridge G A Kuck S Petermann K Pohlmann U Huber G Nozaki H Itoh S 1996 J. Phys. Chem. Solids 57 41
[5] Zhao Y He D W Daemen L L Shen T D Schwarz R B Zhu Y Bish D L Huang J Zhang J Shen G Qian J Zerda T W 2002 J. Mater. Res. 17 3139
[6] Fan X F Wu H Y Shen Z X Kuo J L 2009 Diamond & Relat. Mater. 18 1278
[7] Luo X G Zhou X F Liu Z Y He J L Xu B Yu D L Wang H T Tian Y J 2008 J. Phys. Chem. C 112 9516
[8] Luo X G Guo X J Xu B Wu Q Hu Q Liu Z He J Yu D Tian Y Wang H T 2007 Phys. Rev. B Condens. Matter 76 094103
[9] Chen S Y Gong X G Wei S 2007 Phys. Rev. Lett. 98 015502
[10] Qu N R Wang H C Li Q Li Z P Gao F M 2019 Chin. Phys. Lett. 36 036201
[11] Hou L Gao F M Gou H Y Wang Z B Tian M 2008 Cryst. Growth & Des. 8 1972
[12] Zang C P Sun H 2010 Phys. Rev. B 81 012106
[13] Chen Z Q Li C M Wang J Wang F Wang Q 2013 Sci. Sin. 43 142
[14] Wang D Shi R Gan L H 2017 Chem. Phys. Lett. 669 80
[15] Tang M J He D W Wang W D Wang H K Xu C Li F J Guan J W 2012 Scr. Mater. 66 781
[16] Zhang Y Sun H Chen C F 2008 Phys. Rev. B 77 094120
[17] Li F Man Y H Li C M Wang J P Chen Z Q 2015 Comput. Mater. Sci. 102 327
[18] Chen S Y Gong X G Wei S H 2008 Phys. Rev. B 77 014113
[19] Kawaguchi Wakukawa Kawano 2001 Synth. Met. 125 259
[20] Luo X G Guo X J Liu Z Y He J L Yu D L Tian Y J Wang H T 2007 J. Appl. Phys. 101 083505
[21] Zhang X X Wang Y C Lv J Zhu C Y Li Q Zhang M Li Q Ma Y M 2013 J. Chem. Phys. 138 114101
[22] Ullah S Wang L Li J X Li R H Chen X Q 2019 Chin. Phys. B 28 077105
[23] Wang Y C Lv J Zhu L Ma Y M 2010 Phys. Rev. B 82 094116
[24] Wang Y Lv J Zhu L Ma Y 2012 Comput. Phys. Commun. 183 2063
[25] Li Q Ma Y M Oganov A R Wang H B Wang H Xu Y Cui T Mao H K Zou G T 2009 Phys. Rev. Lett. 102 175506
[26] Lv J Wang Y C Zhu L Ma Y M 2011 Phys. Rev. Lett. 106 015503
[27] Li Q Zhou D Zheng W T Ma Y M Chen C F 2013 Phys. Rev. Lett. 110 136403
[28] Li Y Hao J Liu H Y Lu S Y Tse J S 2015 Phys. Rev. Lett. 115 105502
[29] Hohenberg P Kohn W 1964 Phys. Rev. 136 B864
[30] Kohn W Sham L J 1965 Phys. Rev. 140 A1133
[31] Blöchl P E 1994 Phys. Rev. B 50 17953
[32] Perdew J P Burke K Ernzerhof M 1996 Phys. Rev. Lett. 77 3865
[33] Ding Y 2012 Physica B 407 2282
[34] Lazzeri M Mauri F 2003 Phys. Rev. Lett. 90 036401
[35] N Troullier J L M 1991 Phys. Rev. B 43 1993
[36] Roundy D Krenn C R Cohen M L Morris J W Jr. 1999 Phys. Rev. Lett. 82 2713
[37] Zhang Y Sun H Chen C F 2006 Phys. Rev. B 73 144115
[38] Zhang R F Veprek S Argon A S 2008 Phys. Rev. B 77 172103
[39] Chen X Q Fu C L Podloucky R 2008 Phys. Rev. B 77 064103
[40] Guo W F Wang L S Li Z P Xia M R Gao F M 2015 Chin. Phys. Lett. 32 096201
[41] Gao F M He J L Wu E D Liu S M Yu D L Li D C Zhang S Y Tian Y J 2003 Phys. Rev. Lett. 91 015502
[42] Gao F M Xu R Liu K 2005 Phys. Rev. B 71 052103
[43] Gao F M 2006 Phys. Rev. B 73 132104
[44] Gutiérrez G Menéndez-Proupin E Singh A K 2006 J. Appl. Phys. 99 103504
[45] Wu Z J Zhao E J Xiang H P Hao X F Liu X J Meng J 2007 Phys. Rev. B 76 054115
[46] Pugh S F 1954 Philos. Mag. A 45 823
[47] Yoo C S Kim M Lim J Ryu Y J Batyrev I G 2018 J. Phys. Chem. C 122 13054
[48] Li Z P Gao F M Xu Z M 2012 Comput. Mater. Sci. 62 55
[49] Li Z P Gao F M Xu Z M 2012 Phys. Rev. B 85 144115
[50] Pan Z C Sun H Chen C F 2007 Phys. Rev. Lett. 98 135505
[51] Li Z P Gao F M 2012 Phys. Chem. Chem. Phys. 14 869
[52] Chen X Q Niu H Li D Li Y 2011 Intermetallics 19 1275
[53] Chen X Q Niu H Franchini C Li D Z Li Y Y 2011 Phys. Rev. B 84 121405
[54] Siethoff H 2000 J. Appl. Phys. 87 3301
[55] Phillips J C 1970 Rev. Mod. Phys. 42 317